arXiv:1507.04377vl [cond-mat.soft] 15 Jul 2015 


Submitted, 15.07.2015 


A master equation for force distributions in soft particle packings - Irreversible 
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Mechanical responses of soft particle packings to quasi-static deformations are determined by 
the microscopic restructuring of force-chain networks, where complex non-afUne displacements of 
constituent particles cause the irreversible macroscopic behavior. Recently, we have proposed a 
master equation for the probability distribution functions of contact forces and interparticle gaps 
[K. Saitoh et al., Soft Matter 11, 1253 (2015)], where mutual exchanges of contacts and interparticle 
gaps, i.e. opening and closing contacts, are also involved in the stochastic description with the aid of 
Delaunay triangulations. We describe full details of the master equation and numerically investigate 
irreversible mechanical responses of soft particle packings to cyclic loading. The irreversibility 
observed in molecular dynamics simulations is well reproduced by the master equation if the system 
undergoes quasi-static deformations. We also confirm that the degree of irreversible responses is a 
decreasing function of the area fraction and the number of cycles. 

PACS numbers: 45.70.Cc,46.65.+g,61.43.-j 


I. INTRODUCTION 


Soft particle packings , e.g. colloids, emulsions, foams, 
glasses, and granular materials, are ubiquitous in nature 
and a better understanding of their mechanical proper¬ 
ties is crucial for industry and science [lj. Different from 
usual solids, their constituents are macroscopic particles 
such that thermal fluctuations are negligibly small for 
the individual motions (except for glasses [2J) and thus 
their macroscopic behavior purely originates from the 
mechanics of constituent particles Q. Apart from some 
crystalline systems, e.g. colloidal crystals (dj], their con¬ 
figurations are mostly random (or in amorphous states) 
so that they are temporarily at rest, in mechanical equi¬ 
librium, once the system has relaxed to a static state 
EWul- Then, a packing fraction , <fi , has been used as a 
measure of rigidity of soft particle packings. Some phys¬ 
ical quantities responsible for either mechanical proper¬ 
ties or microscopic structure exhibit the critical behav¬ 
ior near the onset of loss of rigidity (T3 - [Tg[ : The static 
pressure, shear modulus, and excess coordination num¬ 
ber vanish at the onset, while the first peak of the radial 
distribution function diverges near the rigidity transition 
HSUS (such a divergence is specific to static packings, 
which is smoothed out once temperature is imposed to 
the system [2ll - l23| ). In addition, normal mode analy¬ 
ses have revealed excessive low-frequency modes or soft 
modes near the transition, implying large-scale collective 
motions of constituent particles without any changes of 
elastic energy [24l-|27[. Then, introducing the onset as 
a jamming packing fraction , cf>j , a wide variety of ther¬ 
mal (e.g. glasses) and athermal (e.g. granular materials) 
soft particles is unified in a phase diagram [28l - [30j] , while 
there are still some discussions about finite size effects 
[3l| , f32| and the uncertainty of cfj [33 . 

Mechanical properties of soft particle packings have 
their microscopic origin in complex networks of interpar¬ 


ticle (or contact) forces, i.e. force-chain networks [341 ]. 
Therefore, any macroscopic quantity can be deduced 
from statistical averages over the probability distribu¬ 
tion functions (PDFs) of contact forces, which have 
been widely investigated through experiments 
and molecular dynamics (MD) simulations of friction¬ 
less mM 01 frictional particles [48j - [5l[ . In general, 
the PDFs are asymmetric and cannot be described by 
conventional distribution functions [52| such that there 
is still much debate about their tails [53j and asymp¬ 
totic behavior near the zero contact force [53 - l57| . As 
a result, there have been many attempts to establish a 
statistical mechanics of “static” particle packings. Here, 
the main idea is that static configurations for the same 
packing fraction are assumed to be equiprobable, where 
the density of states or PDF of forces is to be calculated 
from appropriate ensembles satisf ying some mechanical 
constraint s, e . g. E dwards’ entropy 
imization _[66h72| , force ensembles 
bles [79M8ll ] . or others [821483 ] . 

When global deformations or affine deformations [85l] 
are applied to soft particle packings, however, the parti¬ 
cles rearrange to rest on other stable (and more favorable) 
configurations. Such rearrangements, i.e. non-affine dis¬ 
placements :, of constituent particles, cause anomalous 
mechanical responses of soft particle packings, espe¬ 
cially, irreversible responses to quasi-static deformations 
lj, Ei2l IHHII1- The degree of non-affinity tends to be 
strong near the jamming transition [8914911 ] and vortex¬ 
like structures of non-affine displacements have been ex¬ 
tensively invest igate d by experiments [p^ iiMl ] and MD 
simu lations with the focus on anisotropy fioTl - 

Il03l ]. As expected, spatial distributions of large non- 
affine displacements can be ma pped onto localized struc¬ 
tures of low-frequency modes 104l - ll07[ . which induces 
notable ano malie s of local clastic constants of soft parti¬ 
cle packings [I08l4ll0l ]. 
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During such non-affine deformations, we observe the 
complicated restructuring of force-chain networks involv¬ 
ing recombinations of contacts, where existing contacts 
are broken and new contacts are generated, i.e. opening 
and closing contacts , respectively. Then, non-affine re¬ 
sponses of contact forces can be regarded as stochastic 
(rather than deterministic affine responses, e.g. effective 
medium theory [85}). Therefore, it is natural to describe 
the evolution of the PDFs by a stochastic model as a 
consequence of stochastic changes of force-chain networks 
during non-affine deformations. Recently, we have pro¬ 
posed a master equation for the PDFs [881] . where the 
master equation fully describes the non-affine evolution 
of the PDFs and anomalous responses of macroscopic 
quantities during deformations. 

In this paper, we explain full details of the master equa¬ 
tion and investigate further applications with the focus 
on irreversible responses of soft particle packings. In the 
following, we introduce our numerical model in Sec. [TT] 
and show our results in Sec. EH Then, we discuss and 
conclude our results in Sec. m and give some technical 
details in Appendices lAl and iBl 


II. METHOD 


We use molecular dynamics (MD) simulations of two- 
dimensional 50 : 50 binary mixtures of frictionless soft 
particles with the same mass, m, and two kinds of radii, 
Ri and Rj ( Ri/Rj = 1.4). The normal force between the 
particles in contact is given by fij = kxij + rjiij with a 
spring constant, k, and viscosity coefficient, rj. Here, Xij 
is an overlap between the particles defined as 


X i — Ri 


Rj dij 


(1) 


with an interparticle distance, dij , and Xij is the relative 
speed in the normal direction. A global damping force, 
f'lamP = —r]Vi, proportional to the particle’s velocity, 
Vi, is also introduced to enhance the relaxation, where 
the damping coefficient is the same with the viscosity 
coefficient between the particles in contact. 

To make static packings of N particles, we randomly 
distribute them in a L x L square periodic box, where no 
particle touches others. Then, we rescale every radius as 


Ri{t + St) 


! x - x m (t) 


Ri(t) 


( 2 ) 


if the averaged overlap is smaller (larger) than the tar¬ 
get value, x m (t) < x (x m (t) > x), so that the averaged 
overlap will finally converge to x. Note that the ratio be¬ 
tween different radii does not change by the rescaling, i.e. 
Ri (t + St) / Rj (t + St) = Ri (t) / Rj (t ). We then stop the 
rescaling when every acceleration of particles drops be¬ 
low a threshold, 10~ 6 kd/m, and assume that the system 
is static. 

We apply an isotropic compression or decompression 
to the prepared packings by rescaling every radius as 


R'i = \j 1±S j R i ( 3 ) 

so that the area fraction increases or decreases from cj> to 
cj)±S(j). Then, we relax the system until every acceleration 
of particles drops below the threshold again. 

In our simulations, distances from jamming are deter¬ 
mined by the linear scaling of averaged overlap, x(4>) = 
A (<t> - <t>j) El, where we estimate the jamming den¬ 
sity as <j>j = 0.8458 ± 10~ 4 with the critical amplitude, 
A = (0.9 ± 0.003)ct, from our 50-sample simulations of 
N = 8192 particles. We also prepared 20 samples for 
smaller systems (N = 512, 2048) and 2 samples for the 
largest one (N = 32768) by changing the initial configu¬ 
rations, while we mainly report the results of N = 8192 
since no result depends on the system size (see Fig. E3- 


III. RESULTS 

In this section, we introduce a master equation for 
the PDFs of particle overlaps. The master equation de¬ 
scribes microscopic changes of force-chain networks dur¬ 
ing quasi-static deformations and can be connected to 
the constitutive relations through the derivative of the 
PDFs. We first study microscopic responses of force- 
chain networks to quasi-static isotropic (de) compressions 
ISec. lIII All and then introduce a master equation for the 
PDFs (Sec. IIII Bl) . To formulate the master equation, we 
quantify mean values and fluctuations of overlaps (Sec. 
IIII Cl) and numerically determine transition rates in the 
master equation (Sec. IIII Dl) . We validate our framework 
by comparing numerical solutions of the master equation 
with MD simulations (Sec. IIII El) , where irreversible re¬ 
sponses of soft particle packings to cyclic loading are also 
examined (Sec. IIIIFl ). 


(i = 1 ,...,N), where t, St, x, and x m (t) are time, an 
increment of time, a target value of averaged overlap, 
and the averaged overlap at time t, respectively. Here, 
we use a long length scale l = 10 2 <r with the mean di¬ 
ameter, er, to rescale each radius gently. We confirmed 
that static packings prepared with longer length scales, 
l = 10 3 tr and 10 4 ct, give the same results concerning crit¬ 
ical scaling near jamming , while we cannot obtain 
the same results with a shorter length scale, l = 10 a. 
During the rescaling, each radius increases (decreases) 


A. Microscopic responses 

At microscopic scales in soft particle packings, mechan¬ 
ical responses to quasi-static deformations are probed as 
restructuring of force-chain networks, where complicated 
particle rearrangements cause the recombination of force- 
chains, i.e. opening and closing contacts. To take into 
account such opening and closing contacts, we employ 
the Delaunay triangulation (DT) of particle packings as 
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FIG. 1. (Color online) The Delaunay triangulation (DT) of 
a soft particle packing, where the red and blue solid lines 
connect the particles in contacts and virtual contacts, respec¬ 
tively. The width of red solid lines is proportional to the 
strength of the interparticle force, where the number of par¬ 
ticles is N = 512. 


FIG. 2. (Color online) Non-affine displacements of IV = 8192 
particles (the red arrows), where the color coordinate (from 
0 to 1) represents the magnitude of non-affine displacements 
scaled by the maximum value. 


shown in Fig. [TJ where not only the particles in contacts, 
but also the nearest neighbors without contacts, i.e. the 
particles in virtual contacts , are connected by the Delau¬ 
nay edges. We then generalize overlaps, Eq. m , as 

Xij = Ri “t“ Rj Rij (d) 

with the Delaunay edge length, D,;j . where the overlaps 
(or gaps) between particles in virtual contacts (Ri+Rj < 
Dij) are defined as negative values. Because the DT is 
unique for each packing, virtual contacts are uniquely 
determined. 

If we apply isotropic compression or decompression to 
the system by Eq. ([3]), every generalized overlaps (not 
only contacts, but also virtual contacts) changes to 

(5) 

where we neglect ed th e higher order terms proportional 
to XijStfi and Sc /) 2 [ml. However, particles are randomly 
arranged and each force balance is broken by the affine 
deformation. Then, the particles move and the system 
relaxes to a new equilibrium state, where non-affine dis¬ 
placements during the relaxation (Fig. [2]) cause compli¬ 
cated ch anges of contacts, including opening and closing 
contacts jl!2| . 

After the relaxation, overlaps change to new values, 
x'ij 7 ^ x i^ ne , that is non-affine responses of overlaps. As 
shown in Fig. 01 there are only four kinds of changes 
from Xij to x[j\ A positive overlap, Xij > 0, remains pos¬ 
itive, xtj > 0, or a negative overlap, < 0, stays neg¬ 
ative, x\j < 0 , such that contacts are neither generated 
nor broken. We call these changes “ contact-to-contact 
(CC)” and “ virtual-to-virtual (VV)” transitions, respec¬ 
tively. On the other hand, if a positive overlap changes to 
a negative one and a negative overlap becomes positive, 
an existing contact is broken and a new contact is gener¬ 
ated, respectively. We name these changes “ contact-to- 
virtual (CV)” and “virtual-to-contact (VC)” transitions, 
respectively. 

In the following, we scale the generalized overlaps by 
the averaged overlap before deformation, x((f>). Then, the 
affine response, Eq. m , is scaled as 

£ affine = £ ± Bal ( 6 ) 


(a) 



(l) ) (VC) (CC) 



FIG. 3. (Color online) Sketches of contacts (a) before and 
(b) after deformation, where four kinds of transitions are 
displayed: (CC) contact-to-contact, (VV) virtual-to-virtual, 
(CV) contact-to-virtual, and (VC) virtual-to-contact, respec¬ 
tively. 


which is a linear function of a scaled overlap before de¬ 
formation, £ = Xij/x(<f>). Here, we omit the subscript, ij, 
after the scaling. On the right-hand-side of Eq. ©, the 
offset is proportional to a scaled strain increment , 


8<t> 

1 ~ ’ 


(7) 


with the amplitude defined as B a = /(2A<j>). Simi¬ 

larly, the non-affine response is scaled as £' = x'^/x^) ^ 

^affine 

Note that we only analyze contact changes occurring 
on the generalized force-chain networks. If the applied 
strain increment is too large, we also observe that parti¬ 
cles, which were neither in contact nor in virtual contact, 
are connected by Delaunay edges after deformation, and 
vice versa. However, we confirm that such rare events 
are below our statistical significance for the whole range 
of sc aled strain increments used in our MD simulations 
m- As shown in Fig. 01 the probability of finding vir¬ 
tual contacts, which are broken to or generated from the 
pairs neither in contacts nor in virtual contacts, i.e. £vn 
or £nv, is less than 3% (the open squares and circles). 
Similarly, the probability of finding contacts which are 
broken to or generated from such pairs, i.e. £qn or £nc, 
is less than 0.1% (the open triangles). Though the two 
cases, (CC) and (VV), dominate the number of contact 
changes, the number of closing and opening contacts, 
(VC) and (CV), asymptotically decays to zero with the 
scaled strain increment as £vc,£cv ~ 7°' 7 if 7 < 1 . 
Therefore, the system exhibits non-affine deformations 
even if 7 is very small. 


B. A master equation 

The restructuring of force-chain networks attributed to 
the contact changes, (CC), (VV), (CV), and (VC), is well 
captured by the PDFs of scaled overlaps, P^(£), where 
the subscript represents the area fraction. Because the 
total number of contacts and virtual contacts is conserved 
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FIG. 4. (Color online) Probabilities of contact changes, £ a p, 
plotted against the scaled strain increment, 7 , where we take 
50-sample averages of IV = 8192 particles for each value of 7 , 
i.e. for each combination of 5(p and <j> — (pj. The subscripts, 
a and /3 (= C,V,N), represent the status before (a) and 
after (/3) compressions, where C, V, and N mean that the 
particles are in contact, in virtual contact, and in neither 
contact nor virtual contact, respectively. Different symbols 
represent different contact changes as listed in the legend. 
The two dashed lines are power law fits for the probabilities, 
£vc - 3.17 x 10 ~ 2 7 0 - 7 and £ C v - 0.94 x HrV 7 . 



-6 -3 0 3 6 

scaled overlap 


FIG. 5. (Color online) The PDFs of scaled overlaps, P ,),(£) 
(squares), P^+s^ (Affine) (triangles), and P^+s^') (circles), 
obtained from our MD simulations of N = 8192 particles 
with the distance from jamming, <p — <pj = 1.2 x 10~ 3 . Here, 
we apply an isotropic compression with S(p = 4 x 10” 4 , i.e. 
7 = 5<f>/ (cp — (pj) = 0.33. The inset is the zoom-in to the 
PDFs of virtual contacts, where the error-bars are obtained 
from 50-sample averages. Because of the bidispersed diam¬ 
eters of particles, the discontinuous gap around zero in the 
initial PDF, P^(£), is smoothed out after the affine deforma¬ 
tion, P(/)-\-8(l> (£affine) - 


during deformations, the PDFs are normalized as 

/ OO 

mn =1 • (8) 

-00 

Figure [5] shows the PDFs obtained from our MD sim¬ 
ulations before compression, P^(£), after affine defor¬ 
mation, P 0 + ^ 0 (^ affine ), and after non-affine deformation, 
P<j>+s<t>{£,')i where the affine compression just shifts the 
initial PDF to the positive direction, while the non-affine 
deformation broadens it in positive overlaps and main¬ 
tain a discontinuous “jump” around zero [22|. Note that, 
however, the PDF after non-affine deformation in nega¬ 
tive overlaps is comparable to that after affine deforma¬ 
tion (see the inset in Fig.0. 

To describe such the non-affine evolution of the PDFs, 
we assume that transitions between overlaps (from £ to 
£') are Markov processes such that we can connect the 
PDF after non-affine deformation to that before compres¬ 
sion through the Chapman-Kolmogorov equation 11^, 

/ OO 

wg . ( 9 ) 

-OO 

where 1 / F(^ , |^) is the conditional probability distribution 
(CPD) of the overlaps, £', which were £ before the com¬ 
pression. In the Chapman-Kolmogorov equation ©, the 
source- and sink-terms caused by the rare events, £yi v, 
£nv ■, £cn, and £nc , are not taken into account. By 
definition, the CPD is normalized as |l!4j 

nOO 

/ wwod? = 1 . ( 10 ) 


From Eqs. ([9l) and (HOD , a master equation for the PDFs 
is readily found to be |l!4j ] 

= I- ^'^ p ^- T ^ p ^ d ^ > 

(11) 

if we introduce a transition rate |114j | as 

T(£'|£) = lim• ( 12 ) 

<50->O 0(p 

The first and second terms in the integral on the right- 
hand-side of the master equation m represent the gain 
and loss of overlaps, £', respectively. Therefore, the tran¬ 
sition rates or CPDs fully determine the statistics of mi¬ 
croscopic changes of force-chain networks. 

If we multiply Eq. CEB by (p — (pj and introduce an 
infinitesimal scaled strain increment as dy = 5<p/((p — 
cpj ) < 1 , we obtain an alternative form of the master 
equation as 

|^(0 = [r 7 (£'|£)P*(£) - T 7 (£|£')P*(£')] , 

(13) 

where the transition rates are now defined as P 7 (£ , |£) = 
lim^o W(£'|£)/<57- 

In addition, the master equation for unsealed overlaps, 
x = x {</>)£ and x' = x(cp)^', can be written as 

8 r°° 

^V) = J_ [T*(x'\x)p;(x) - T*(x\x')p;(x')\ dx , 

(14) 


— OO 
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where the unsealed PDF, unsealed transition rate, and 
unsealed increment are given by P£(x) = i(0) _1 P0(£), 
T*( x'\x) = i(</>)- 1 T(£ , |£), dx = x((p)d£, respectively. 

The PDFs can be connected to macroscopic quantities 
through the n-th moment of positive overlaps, 

pOO pOO 

p n = / x n P <j> {x)dx = / x n P (t> {£ > )df , (15) 

Jo Jo 


where we used the relation, P l p(x)dx = P(/,(f)df. Then, 
the coordination number, z, averaged overlap, x, and 
static pressure, p, are given by the first three moments 
as 


2 N e 

z = ~l \T Mo ’ 
x = Mi , 
kN E 

p = wmi - P2> , 


(16) 

(17) 

(18) 




£ f 


respectively, where iV# is the total number of the Delau¬ 
nay edges [l 15j | . 


C. Mean and fluctuations 

Because the transition rates or CPDs are defined as 
distributions of Ef around their mean values, we first mea¬ 
sure the mean and fluctuations of £' through scatter plots 
of scaled overlaps. Figure [ 6 ] displays the scatter plots ob¬ 
tained from our MD simulations under compressions. In 
this figure, the four kinds of transitions displayed in Fig. 
[3] are mapped onto four regions: (CC) £,£' > 0, (VV) 
< 0, (CV) £ > 0, £ < 0, and (VC) f < 0, ? > 0, re¬ 
spectively. Though scaled overlaps after affine deforma¬ 
tion, £ affine , are described by the deterministic equation 
ED, those after non-affine deformation distribute around 
their mean values with finite fluctuations, i.e. non-affine 
responses of overlaps are stochastic. The differences be¬ 
tween affine and non-affine responses are always present, 
but not visible if the applied strain is small or the sys¬ 
tem is far from jamming, i.e. if 7 <C 1 (Fig. (HDa)), while 
deviates more from £ affine and data points are more 
dispersed if 7 1 (Fig. (Hid)). 

In (CC) and (VV) regions, the mean values of £' can 
be fitted by linear functions of 

m(Q = i a i + 1)£ + b i , (19) 

where the subscripts, l = c and v, represent the mean 
values in (CC) and (VV), respectively. We also introduce 
standard deviations of £' from their mean values as Vi, 
which are found to be almost independent of £. Then, the 
systematic deviation from the affine response, £ affine , is 
quantified by the coefficients, a/, 5/, and vi as summarized 
in Fig. [ 8 ] Note that the affine response, Eq. ED, is 
recovered if ai = vi = 0 and bi = B a 7 . 

As we observed in the scatter plots (Fig. [G]), the dif¬ 
ference between affine and non-affine responses increases 


FIG. 6 . (Color online) Scatter plots of scaled overlaps under 
compression, where the blue and red dots are £ affine and (f 
plotted against £, respectively. Distances from jamming are 
<j> - (j>j = (a) 4 x 10" 3 , (b) 1.2 x 1CT 4 , (c) 1.2 x 1(T 3 , and 
(d) 4x 10 ~ 4 , respectively, while applied strain increments are 
8(j> = 4 x 10 -5 ((a) and (b)) and 4 x 10 -4 ((c) and (d)), 
respectively, such that scaled strain increments are given by 
7 = (a) 0.01, (b) 0.33, (c) 0.33, and (d) 1, respectively. Here, 
the number of particles is TV = 8192 (1 sample). 




£ £ 


FIG. 7. (Color online) Scatter plots of scaled overlaps under 
decompression, where the blue and green dots are £ affine and 
plotted against £, respectively. Applied strain increments 
are <50 = 4 x 10 ~ 5 ((a) and (b)) and 4 x 10 -4 ((c) and (d)), 
respectively, where distances from jamming, magnitudes of 
scaled strain increments, and the number of particles are as 
in Fig. [0 (d) If 7 > 1, an unjamming transition happens. 
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FIG. 8 . (Color online) A schematic picture of the difference 
between affine and non-affine responses of scaled overlaps, 
where the blue and red solid lines represent £ affine (for small 
and large particles) and linear functions, mi(t ;), respectively. 
In (CC), £ affine and rrii(£) intersect at £ = £*, where the green 
arrows represent the decrease and increase of scaled overlaps 
during relaxation. The excess slope in (CC), a c , and all the 
dimensionless lengths, bi, vi, and A i, are proportional to 7 , 
where A c and A„ represent typical penetration lengths of new 
contacts and new virtual contacts, respectively. 

with the scaled strain increment, 7 . From our simula¬ 
tions, we find that all the coefficients, except for a v — 0 , 
linearly increase with 7 , where all data with a wide va¬ 
riety of 8(f) and <f> — (f>j collapse onto a linear scaling of 7 
(Fig.[U)a)). As shown in Figs. [91(b)- (f), all the coefficients 
(the open squares) are well described by linear scalings 
of the scaled strain increment, 

ai = Aij , b t =Bi"/, vi=Vi\j\, (20) 

where the scaling amplitudes, Ai, Bi , and V), estimated 
in the fitting range, ICC 6 < 7 < 1, are listed in Table 
[H Because a v ~ 0 and B v ss -B a (— 1.3 for small and 
large particles), virtual contacts almost behave affine in 
average except for their huge fluctuations ( V v V c ). 

In contrast, B c is always smaller than B a such that 
m c (£) intersects £ affine at a characteristic scaled overlap, 
£* = ( B a — B c )/A c ~ 1.4 (which is independent of 7 ). 
This leads to small responses (£' < £ a ffine) of small over¬ 
laps (£ < £*) and vice versa as indicated by the green 
arrows in Fig. [51 implying preferred tangential and hin¬ 
dered normal displacements as a sign of non-affine defor¬ 
mations [90]. As shown in Fig. [TUI the linear scalings, 
Eq. (1501) . do not depend on the system size. 

The linear scalings, Eq. 0 , are retained under de¬ 
compression, i.e. for — 7 . Figures [7] displays scatter plots 
of scaled overlaps under decompression, wher e an unjam¬ 
ming transition happens if 7 > 1 (Fig.[7](d)) jlld |. Here, 
we also fit the mean values of £' by m;(£) and quantify 
the fluctuations by standard deviations, 17 . Then, we find 
that the excess slope and offsets are negative, a Cl bi < 0 , 
while a v ~ 0 and vi > 0 as in the case of compression. As 







FIG. 9. (Color online) (a) Data collapse of the excess slope in 
(CC), a c , plotted against 7 = 5(f>/{(f> — (f>j) under compression, 
where different symbols represent 8(f) as listed in the legend. 

(b) -(f): Double logarithmic plots of the coefficients, (b) |o c |, 

(c) | b c I, and (d) | 6 „|, and standard deviations, (e) v c and (f) 
v v , under compression (squares) and decompression (circles), 
where we averaged the data over different combinations of 5(f) 
and (f> — (f>j for the same 7 . Here, the dashed lines represent 
the linear scaling against 7 , i.e. Eq. (1201) (the vertical dotted 
lines represent error-bars comparable to the mean values). 

shown in Figs. [HKb)-(f), their absolute values (the open 
circles) are well described by the same linear scalings for 
compression so that Eq. (l20l) can be used for both com¬ 
pression and decompression (because a c and bi are linear 
against 7 , i.e. a c ,bi ~ 7 , they are positive and negative 
under compression, 7 , and decompression, — 7 , respec¬ 
tively). 

In contrast to (CC) and (VV), the data of £' in (VC) 
and (CV) are concentrated in narrow regions (the inside 
of the dashed lines in Fig. [51 left and below, respectively). 
Here, £ affine linearly increases with £ in (VC) and there 
is no data of £ affine in (CV), because affine compressions 
do not generate any opening contacts. Similarly, affine 
decompressions do not generate any closing contacts. 

D. Conditional probability distributions 

The statistics of restructuring of force-chain networks 
is fully described by the transition rates in the mas- 
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TABLE I. Scaling amplitudes in Eqs. ((201) . q-indices in Eq. (Ei, and dimensionless length scales in Eqs. (f27|) and (12811 . 


l 

A, 

B l 

Vi 

Qi 

Az 

c 

0.76 ±2.4 x lCT 11 

0.24 ±3.3 x 10“ 3 

0.32 ±4.8 x 10 -3 

1.13 ±0.13 

6.10 ± 0.33 

V 

0.00 ±8.3 x 10“ 3 

1.80 ±8.4 x 10 -2 

4.41 ± 1.9 x 10 _1 

1.39 ±0.18 

0.65 ± 0.04 
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FIG. 10. (Color online) The system size dependence of the 
coefficients under compression, (a) |a c |, (b) a„, (c) | 6 C |, (d) 
| 6 „|, (e) v c , and (f) v v , where the dotted lines represent the 
linear scalings of 7 , i.e. Eq. (ESI, except for a v ~ 0. The 
different symbols represent different system sizes, N, as listed 
in the legend in (a). 


ter equation m or CPDs in the Chapman-Kolmogorov 
equation ©■ For example, the CPD for affine deforma¬ 
tion is given by a delta function as 

VFaffi„e(e'IC)=^'-r ffine ) , (21) 

which just shifts the PDF by -B a 7 , be. P^+s^iO = 
Prj ,(£ — Ba'y), as shown in Fig. [5] However, the CPDs 
for non-affine deformations distribute around the mean 
values, mj(£), with finite widths, vi. In the following, we 
determine the CPDs after non-affine deformations from 
our results of MD simulations, where we explain how to 
calculate the CPDs from numerical data in Appendix [A] 
First, we determine the CPDs in (CC) and (VV) un¬ 
der compression. Figures [III a) and (b) show the CPDs 
in (CC) and (VV), respectively, where all data are sym¬ 
metric around the mean values and well collapse if we 


multiply the CPDs and distances from the mean values, 
Ei = £' — mi(£), by 7 and 1 / 7 , respectively. In these 
figures, the CPDs in (CC) and (VV) are well described 
by the solid lines, 


7Wcc(£'|£) = /c(S c /7) , (22) 

iWw^M) = /«(S„/7) , (23) 

respectively, where 


fi(x) 


1 

c(qi) 


1 + 



1 -q t 


(24) 


is the q-Gaussian distribution [U7tini . InEq. the 
functions are introduced as n(t) = (t — 3) / (1 — t) and 
c(t) = Vi\fn{t)B (1/2, n{t)/2) with the beta function, 
B(x,y). In Tabic HI we list the q-indices (1 < qi < 3) 
which determine shapes of the CPDs, where the normal 
(Gau ssian) distribution is recovered in the limit of qi —> 1 

Next, we determine the CPDs in (CV) and (VC) un¬ 
der compression, where the normalization condition, Eq. 
m, is satisfied in £ > 0 and £ < 0 as 


poo p0 

/ W cc {t'\€)d?+ WcviZ' 10^ = 1, (25) 

J 0 J —00 

poo p 0 

/ w vc (?m'+ / WW(£'|£K' = 1, (26) 

J 0 J —00 


respectively. Figures nn c) and (d) show the CPDs in 
(CV) and (VC), respectively, where all results collapse 
after the same scaling as for those in (CC) and (VV). 
In these figures, the CPDs in (CV) and (VC) are well 
described by the exponential distributions , 


7W / cf(£ / |£) = Fc 


ra c (g) 

v c 


p«'/7A„ 


Xv 


jWvc(,€' I?) = Fy 


m v {£) 

v v 


e ~P/j\c 

A c 


(27) 

(28) 


respectively, where the dimensionless length scales, A j, 
are listed in Table Q] On the right-hand-sides of Eqs. (EH 
and ESI , the ^-dependent functions are defined as the 
incomplete beta function, 


Fi(x) = -B 


i (qi) 


Mi) 


n{qi) 1" 

2 ’ 2 ’ 


(29) 


which are equivalent to the cumulative distributions of 

















































the CPDs in (CC) and (VV), i.e. 


Fc 

’"l c (0 

POO 

= 1 - / WccCe'IO^' , 

(30) 


L v c J 

Jo 


F v 

1 1 

Uy> 

?> 

= 1 - [ Wwieme , 

J —00 

(31) 


respectively. From Eqs. 



e - £ '/TAc 

A c 


dC 


and m, and a relation 


f° g 


I — OO Aj; 


d? = 1 , (32) 


it is confirmed that the distributions defined as Eqs. 
m, nisi), m, and (l28l) satisfy the normalization con¬ 
ditions of the CPDs, Eqs. (E5D and (ESI). As shown in 
Figs. fl2fal and (b), Eqs. (l27l) and (l28l) well describe 
the ^-dependence of the CPDs in (CV) and (VC), re¬ 
spectively. Because the CPDs in (VC) and (CV) ex¬ 
ponentially decay along the £'-axis, the dimensionless 
length scales, A c and A„ (A c -C A„), represent typical 
penetration lengths for new contacts and generated gaps 
between new virtual contacts. Therefore, closing and 
opening contacts contribute to the new PDF, P ( p + s < j ,(£'), 
between — A„ < < A c as shown in Fig. [ 8 ] Note 

that, if 7 = 0, the CPDs defined as Eqs. (1441) . (1451) . 
(1471 ) . and (145D converge to Wcc = Wvv = <5(£ — £') 
and Wcv = Wye = 0 , respectively [l 2 l| . so that the 
Chapman-Kolmogorov equation © does not change the 
PDF without deformation. 

The CPDs under decompression are given by replacing 
the scaled strain increment, 7 , with — 7 , as in the case 
of linear scaling, Eq. for decompression. Figures 

ITTT eWh) display the CPDs under decompression, where 
we multiply the CPDs and distances from the mean val¬ 
ues by 7 and I/ 7 , respectively. Figures fl2T cl and (d) 
show the ^-dependence of the CPDs in (CV) and (VC) 
under decompression, respectively. In these figures, the 
solid lines correspond to Eqs. (1441) . (1451) . (1471) . and (1451) . 
such that the CPDs do not change their shapes under 
compression and decompression. 


E. Numerical validations of the master equation 


We now unfold the master equation for the PDFs of 
scaled overlaps, where the transition rates are divided 
into the four cases, (CC), (VV), (CV), and (VC), as 
Tcc(£ |0 = lim^o WcciZIQ/W, etc. Then, the mas¬ 
ter equation ED can be rewritten for positive and nega¬ 


tive scaled overlaps, > 0 and < 0, as 

= l°° [Tccimp*® - Tccmwo} da 

+ f [Tvciemio-Tcvmwndc, 

J —00 

(33) 

-HjW) = [_ {Tw(e\mo - Twmwo] # 

/»oo 

+ / [Tcv (?moo - Tvcimwe)} , 

Jo 

(34) 

respectively. The second terms on the right-hand-sides 
of Eqs. (1551) and (15(41) are cross-terms of positive and neg¬ 
ative overlaps due to closing and opening contacts. The 
other form of the master equation (1131) can be rewritten 
as well. In Appendix [Bj we derive the master equation, 
Eqs. (1531) and (1541) . from the Chapman-Kolmogorov equa¬ 
tion ©• 

In our framework, transitions between overlaps are as¬ 
sumed to be Markov processes. To validate this assump¬ 
tion, we compare numerical solutions of the master equa¬ 
tion with the PDFs obtained from MD simulations. Fig¬ 
ure El displays the numerical solutions under compres¬ 
sion, where the initial PDF is given by the MD simulation 
with the distance from jamming, <po — cf>j = 4 x 10~ 3 . In 
this figure, the overlaps are scaled by the averaged over¬ 
lap at the initial state, x(tp 0 ), where the increment of area 
fraction is fixed to 8<f> = 4 x 10 ~ 5 such that the scaled 
strain increment is less than 7 = 6 70 = S(f>/((j)o — <j)j) = 
10 -2 . We confirm a good agreement between the nu¬ 
merical solutions (red solid lines) and MD simulations 
(open symbols) even in the tails of the PDFs (the in¬ 
sets in Fig. ED, where the sequential solutions between 
4 x 10 -3 < <j) — <f>j < 4 x 10 -2 are also displayed (dashed 
lines). 

Note that the master equation becomes insensitive to 
the initial condition, which is the only input in our frame¬ 
work, after finite strain steps. Figure El shows numerical 
solutions under compression with different initial con¬ 
ditions, i.e. a step function and Gaussian distribution, 
where the solutions for both initial conditions converge 
to the same PDF after some increments of area fraction. 

The master equation generates (or maintains) discon¬ 
tinuous “jumps” of the PDFs around zero-overlap as 
we have observed after non-affine deformations in Fig. 
[5] Figure El displays numerical solutions of the master 
equation, where the initial PDF is given by a Gaussian 
distribution which is initially continuous around zero. As 
can be seen, a discontinuous jump is generated after sev¬ 
eral increments of area fraction. Figure [151 explains that 
a discontinuous jump is caused by the significant differ¬ 
ence between the dimensionless lengths in the CPDs in 
(CV) and (VC), i.e. A c <C A„: When virtual contacts are 
closed by affine deformation as indicated by the shaded 
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FIG. 11. (Color online) Semi-logarithmic plots of the CPDs in (a) (CC), (b) (VV), (c) (CV), and (d) (VC) under compression, 
and in (e) (CC), (f) (VV), (g) (CV), and (h) (VC) under decompression. Different symbols represent different scaled increments 
as listed in the legends. The solid lines are given by the distribution functions, Eqs. (12211 . (1231) . (1271) . and (1281) . The dotted lines 
in (a) and (e) are Gaussian fits, see Ref. [88j], while those in (b) and (f) are the g-Gaussian distributions for (CC). 


FIG. 12. (Color online) Semi-logarithmic plots of the CPDs 
in (a) (CV) and (b) (VC) under compression, and in (c) (CV) 
and (d) (VC) under decompression. Different symbols repre¬ 
sent different scaled increments as given in Fig. 1111 The solid 
lines are given by the distribution functions, Eqs. (f27|) and 
(1281) , and the insets show the corresponding semi-logarithmic 
plots. 



0 2 4 6 8 10 


scaled overlap 



scaled overlap 


FIG. 13. (Color online) Comparisons between the master 
equation (solid and dashed lines) and MD simulations (open 
symbols) under compression in (a) positive and (b) negative 
overlaps, where the open squares, circles, and triangles are the 
PDFs obtained from our MD simulations at = 4x 10~ 3 , 

1.2 x 10 -2 , and 4 x 10 -2 , respectively. Overlaps are scaled by 
the averaged overlap at (j>o — (f>j = 4 x 10 -3 and the numerical 
solutions develop in the directions indicated by the arrows. 
The insets show the corresponding semilogarithmic plots. 



FIG. 14. (Color online) The dependence of the master equa¬ 
tion on the initial conditions, where the red solid and blue dot¬ 
ted lines are numerical solutions (under compression) starting 
with a step function and Gaussian distribution, respectively. 


FIG. 15. (Color online) The generation of a discontinuous 
jump by the master equation, where the initial PDF (blue 
open squares) is a Gaussian distribution. Numerical solutions 
of the master equation (red open circles) evolve from (a) to 
(d) under compression. 

area in Fig. 1161 a) (no contact is opened by affine com¬ 
pression) , each new contact is pushed back by the repul¬ 
sive interparticle force (the short arrows in Fig. fTCT b) ) 
so that the typical penetration is reduced to A c , which 
is considerably smaller than the maximum penetration 
by affine deformation, B a . Then, the overlaps between 
new contacts decrease and the PDF near zero-overlap 
(on the positive side) grows. However, once particles are 
detached from each other, generated negative overlaps 
can increase freely since there is no interparticle force 
and the typical magnitude of generated negative over¬ 
laps, A v , is much larger than the typical penetration, A c . 
Thus, generated negative overlaps widely distribute such 
that the PDF near zero-overlap on the negative side be¬ 
comes smaller than that on the positive side. As a result, 
a jump between positive and negative sides is generated 
(or maintained). Note that such discontinuities in the 
PDFs are specific to static packings winch will disappear 
once a finite temperature is imposed [HHIIj]. 


F. Irreversible responses 

We next turn to irreversible mechanical responses of 
soft particle packings to quasi-static deformations. 

Figure [T7] shows (a) the coordination number, (b) aver¬ 
aged overlap, and (c) static pressure under compression 
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(a) PDF (b) 



scaled overlap scaled overlap 

FIG. 16. (Color online) Schematic pictures of (a) affine and 
(b) non-affine changes of a PDF. (a) The affine deformation 
shifts the initial Gaussian PDF (thin solid line) to the new 
PDF (blue solid line) as indicated by the arrow, where the 
shaded area corresponds to the amount of new contacts, (b) 
The PDF after affine deformation (blue solid line) changes 
to the PDF after non-affine deformation (red solid line) as 
indicated by arrows, where the typical penetration lengths of 
new contacts and virtual contacts are significantly different, 
Ac ^ . 


PDF 



and decompression, where we first increase the area frac¬ 
tion from (j>Q — <f>j = 4 x 10 -3 to <j>i — (f>j = 8 x 10 -2 
and then decrease it to (f> o — </>j with the increment, 
5(f> = 4 x 10 -4 , assuming that <j)j is not changed by the 
cyclic loading [33}. In this figure, the lines are calcu¬ 
lated by Eqs. (usd-©, where we substituted numerical 
solutions of the master equation for the PDF in the ra¬ 
th moment, Eq. ©. The open symbols are the results 
of MD simulations, where a reasonable agreement with 
the master equation is established. In addition, the mas¬ 
ter equation captures irreversible responses of the macro¬ 
scopic quantities (the coordination number in Fig. I17f a) 
is more visible) and well reproduces the non-linear be¬ 
havior of the static pressure (Fig. ©(c)) 0 - Note that 
the master equation without any opening and closing con- 
tacts , i.e. numerical solutions with zero transition rates 
in (CV) and (VC), gives a linear increase and decrease 
of pressure (straight lines in Figs, ©(b) and (c)) as de¬ 
scribed in some literature focusing on the systems close 
to jamming }l5j j. 


1. Irreversibility of the PDFs 


To quantify the degree of irreversibility, we introduce 
a difference between the PDFs during compression and 
decompression for the same area fraction, P£ fe (£) and 
as 


sk{4>) = 





n 


(35) 


Here, the subscript, k = 1,2,, represents the number 
of cyclic (de)compressions and the superscripts, c and 
d, mean the PDFs during compression and decompres¬ 
sion, respectively. On the right-hand-side of Eq. (l35l) . 
the bracket, (...){, represents the average over all scaled 


overlaps (—oo < £ < oo). Figure ©displays a schematic 
picture of cyclic (de)compressions, where the difference, 
<&(</)), is calculated in the fc-th cycle at the same area 
fraction, <f>. 

First, we study the dependence of irreversibility on the 
scaled strain increment. As shown in Fig. fTTJT a). the dif¬ 
ference in the first cycle, ci (<f>), is a decreasing function of 
the area fraction, where a good agreement between the 
MD simulation (the open triangles) and the master equa¬ 
tion (the blue dotted line) with a small strain increment, 
<5yo = 0.1, is established. Both the differences, 
obtained from MD simulations and the master equation 
increase with the scaled strain increment, where the mas¬ 
ter equation starts to deviate from MD simulations once 
<5yo exceeds unity (Figs, ©(a) and (b)). Figure ©(c) 
shows the dependence of si(<j>) on dyo, where the irre¬ 
versibilities observed in MD simulations (open symbols) 
are well reproduced by the master equation (closed sym¬ 
bols) if the deformation is quasi-static (blue-shaded re¬ 
gion, dyo < 0.1). Remarkably, the difference between 
loading and unloading remains finite, which implies that 
a strong irreversibility present even if the applied strain 
is very small. 

Second, we examine the dependence of irreversibility 
on the number of cyclic (de)compressions. Figure IKiUT a) 
displays the difference, ?/c (</>), for the different number of 
cycles, k, where good agreements between MD simula¬ 
tions (the open symbols) and the master equation (the 
lines) are established with <5yo = 0.1. (In the inset of Fig. 
EUT a). we confirm that the difference decays to zero at the 
turning point, <f> = (f> i.) In this figure, the difference al¬ 
ways decreases with the number of cycles (Figs. l20P a) 
and (b)), where the irreversibilities observed in MD sim¬ 
ulations are well predicted by the master equation (Fig. 
©c)). 


2. Irreversibility and local symmetry 

The microscopic origin of irreversible responses can be 
explained by our observations of the transition rates: If 
the system response is reversible against deformations, 
the transition rates for decompression and compression 
must satisfy a local symmetry, i.e. T_ 7 (£'|£) = T 7 (£|£'). 
This means that the scattered data of scaled overlaps for 
compression and decompression for the same y (Figs. [G] 
and 0 are symmetric with respect to the diagonal line, 

= £. If we assume such a symmetry, however, the mean 
values of £' and standard deviations under decompression 
should be given by m*(£) = £/(aj + 1) — bi/{ai + 1) and 
v* = vi/(ai + 1), respectively, which is clearly different 
from our results, rraj(£) = (1 — ai )£ — bi and vi. Note that 
m*(£) and v* converge to m/(£) and vi in the small strain 
limit, i.e. m*(£) = (1 - a;)£ - bi + 0( y 2 ) ~ m;(£) and 
v* = Vi + 0( y 2 ) ~ vi, so that the scattered data tend to 
be almost symmetric around the diagonal line if y <C 1 . 
On the other hand, the transition rates in (CV) and (VC) 
are finite even when the scaled strain increment is very 
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FIG. 17. (Color online) (a) Coordination number, z, (b) averaged overlap, x, and (c) pressure, p , in units of the mean diameter, 
<7, and spring constant, k, plotted against the area fraction during a compression-decompression cycle (the arrows in (a)). 
The open symbols and lines are obtained from MD simulations and the master equation, respectively. The (red) squares 
and (red) solid lines are the results under compression, while the (blue) circles and (blue) dotted lines are the results under 
decompression. The straight lines in (b) and (c) are given by numerical solutions of the master equation without any opening 
and closing contacts, i.e. Wcv = Wvc = 0 , where the (black) solid and (yellow) dotted lines are the results under compression 
and decompression, respectively. The insets are zooms into the squares surrounded by the broken lines. 


area fraction 



FIG. 18. (Color online) A schematic picture of cyclic 
(de)compressions: During compression, the area fraction in¬ 
creases from 0 o to 0 i along the red arrows, while it decreases 
from 0 i to 0 o along the blue arrows during decompression. 
The difference, 0 ,( 0 ), is defined in the fc-tli cycle at the same 
area fraction, 0, where the PDFs (calculated at the open cir¬ 
cles) are compared with each other. 


small, which leads the finite irreversibility during cyclic 
(de)compressions. In fact, the number of opening and 
closing contacts obeys a power law, Scv^vc ~ 7 °' 7 , as 
observed in Fig. [4j From these observations, the irre¬ 
versible responses cannot be avoided if the applied strain 
is finite (not zero). 


3. Asymptotic behavior near jamming 


We also examine asymptotic behaviors of the PDFs 
near the jamming transition, 0 —> <j)j, by the master 
equation. Here, we numerically solve the master equa¬ 
tion by fixing the scaled increment, where the increment 
of area fraction, 5cj) = (0 — 0 j) 7 , decreases as the sys¬ 
tem approaches to the jamming transition such that the 
system is always above jamming. Then, we introduced a 


difference between the PDFs, P$+s<t>(Q and P^(£), as 

X(0) = , (36) 

where we find that the difference, %( 0 ), asymptotically 
decreases to zero (Fig. [2Tj), implying that the PDF tends 
to be self-similar near jamming. Note that, however, 
the decay is so slow, x(0) ~ (0 — 0j) 014 , and thus it is 
practically impossible to reach the asymptotic limit. 


IV. DISCUSSION AND CONCLUSION 


In summary, we have proposed and numerically inves¬ 
tigated a master equation for the probability distribution 
functions (PDFs) of forces in soft particle packings under 
isotropic compression and decompression. 

First, mechanical responses of soft particle packings 
to quasi-static deformations are determined by the re¬ 
structuring of force-chain networks, involving their com¬ 
plicated recombinations of force-chains, i.e. opening and 
closing contacts. To take into account all kinds of changes 
of contacts, we have introduced Delaunay triangulations 
of soft particle packings, where the force-chain networks 
are generalized to include not only the particles in con¬ 
tact, but also the nearest neighbors in virtual contact. 
Then, the statistics of structuring of force-chain net¬ 
works are well captured by the PDFs of generalized 
overlaps, where the “overlap” (or interparticle gaps) be¬ 
tween the particles in virtual contact is defined as neg¬ 
ative values. In addition, the rates of four kinds of 
changes of contacts, i.e. contact-to-contact (CC), virtual- 
to-virtual (VV), contact-to-virtual (CV), and virtual-to- 
contact (VC), are described by the conditional proba¬ 
bility distributions (CPDs) in the Chapman-Kolmogorov 
equation ©, or equivalently, by the transition rates in 
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FIG. 19. (Color online) (a) The difference for the first cycle, a ((ft), plotted against the area fraction, where the open symbols 
and lines are obtained from MD simulations and the master equation, respectively. Different colors represent different strain 
increments defined at the initial state, i.e. Sjo = <5<(>/(<)>o ~ </>j) = 0.1, 1, and 10, as listed in the legend, (b) The difference, 
a (</>), obtained from the master equation, where the scaled strain increment, < 570 , increases as indicated by the arrow and listed 
in the legend, (c) The dependence of a (4>) on ^7o, where the open and closed symbols are the results of MD simulations and 
the master equation, respectively. Different symbols represent different area fractions as listed in the legend. 



FIG. 20. (Color online) (a) The difference for the k -th cycle, a,(0), plotted against the area fraction, where we use 67,0 = 0.1 and 
the open symbols and lines are obtained from MD simulations and the master equation, respectively. Different colors represent 
the different number of cycles, k, as listed in the legend. The inset shows a zoom-in around the maximum area fraction, <f> 1 . (b) 
The difference, ?*,(<)>), obtained from the master equation with <570 = 0.1, where the number of cycles, k, increases as indicated 
by the arrow and listed in the legend, (c) The dependence of Cfc(<)>) on k, where the open and closed symbols are the results 
of MD simulations and the master equation, respectively. Different symbols represent different area fractions as listed in the 
legend. 


FIG. 21. (Color online) The asymptotic decrease of the dif¬ 
ference between the PDFs, x(0)> plotted against the distance 
from jamming, <j> — (j>j (the open circles). Here, we fixed the 
scaled strain increment to 7 = 0 . 1 , i.e. increments of area frac¬ 
tion are in the range, 10 -12 < 5(f) < 10 -4 . The dotted line 
represents a power law fitting, xifi) = 7.6 x 10 ~ 4 (< j> — (j>j) 01A . 
Note that the initial data drop around <j> — cj>j ~ 10 -3 (the 
open circles in the right-top) is caused by the fluctuations of 
the initial PDF, which is given by the MD simulation. 


the master equation m■ We have numerically deter¬ 
mined the transition rates for the four kinds of changes 
by MD simulations of two-dimensional bidispersed fric¬ 
tionless soft particle packings under both isotropic com¬ 
pression and decompression. 

Second, the transition rates are defined as distribu¬ 
tions of scaled overlaps around their mean values, where 


we have found that the mean values are described by lin¬ 
ear functions of the scaled overlaps before deformation. 
Then, the deviations from affine deformation, character¬ 
ized by the excess slopes (a;) and offsets ( 6 /), linearly in¬ 
crease with the scaled strain increment, 7 = S(f>/(<t> — <j>j), 
for both isotropic compression and decompression. The 
scaled overlaps after non-affine deformation have finite 
widths (vi) around their mean values, where the widths 
also linearly increase with 7 . This is the evidence of 
stochastic processes of overlaps and interparticle gaps 
during the restructuring of force-chain networks. All the 
linear scalings against 7 are not affected by the system 
size. Note that affine deformations do not generate such 
a scatter so that their transition rates are given by delta- 
functions, resulting in a deterministic evolution of the 
PDFs. 

Third, we have found that the transition rates for 
(CC) and (VV) are symmetric around the mean values 
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even though the PDFs are asymmetric and cannot be 
described by conventional distribution functions. The 
transition rates can be unified to the g-Gaussian distri¬ 
butions, where their tails are wider than normal Gaus¬ 
sian distributions implying that the stochastic evolution 
of scaled overlaps is correlated. Note that the tail of the 
transition rate in (VV) is much wider than that in (CC), 
i.e. the correlation between interparticle gaps, interest¬ 
ingly, seems to be much stronger than that between con¬ 
tacts. We have also confirmed that the transition rates 
for (CV) and (VC) are given by exponential distribu¬ 
tions satisfying normalization conditions for £ > 0 and 
£ < 0, respectively. In the transition rates for (CV) and 
(VC), the dimensionless length scales (A i) represent typ¬ 
ical penetration- and gap-lengths after creating contacts 
and virtual contacts, respectively, where their significant 
difference (A c <C A„) induces the discontinuous jump of 
the PDF around zero-overlap. 

Fourth, we have validated the master equation, in 
which stochastics of scaled overlaps are assumed to be 
Markov processes and the initial PDF is the only input, 
by comparing numerical solutions of the master equa¬ 
tion with the PDFs obtained from MD simulations. We 
have confirmed that the evolution of the PDFs is well 
described by the master equation and their dependence 
on the initial condition vanishes with increasing strain 
(deformation). We have also demonstrated the genera¬ 
tion of a discontinuous jump of the PDF by numerically 
solving the master equation. 

Finally, irreversible responses of macroscopic quanti¬ 
ties, e.g. coordination number, averaged overlap, and 
static pressure, defined by the moments of scaled overlaps 
are well reproduced by the master equation. Introducing 
the difference between the PDFs under cyclic compres¬ 
sions and decompressions for the same area fraction, we 
have confirmed that the degree of irreversibility observed 
in MD simulations is well reproduced by the master equa¬ 
tion if the system undergoes quasi-static deformations. 
In both MD simulations and the master equation, we 
have found that the difference is a decreasing function of 
the area fraction and the number of cyclic deformations. 
We have also confirmed the self-similarity of the PDFs by 
the master equation near the jamming transition, where 
the difference between the PDFs asymptotically decays 
to zero as the area fraction approaches to the jamming 
area fraction, (f>j. 

In our recent study [ 1 22t , we have confirmed the basic 
properties of the master equation, i.e. linear scalings for 
the excess slopes, offsets, and fluctuations against 7 , and 
the symmetry of transition rates, by MD simulations of 
two-dimensional polydisperse frictional particles, where 
we have found that the increase of polydispersity and 
decrease of friction broaden the tails of transition rates, 
implying the increase of correlations of contacts and in¬ 
terparticle gaps. In addition, the transition rates have 
been found to be symmetric g-Gaussian dist ribu tions in 
our recent experiments of wooden cylinders [ 1231] . 

Because the master equation requires the increment S<f> 


to be much smaller than (/> — (j)j, i.e. 7 <C 1 , it can never 
reach (f>j. This means that the jamming transition is a 
singular limit of the master equation, but possibly avail¬ 
able to asymptotic analysis. In addition, there is the need 
of further studies on the physical origin of the stochastics 
of overlaps described above. Analytic or asymptotic solu¬ 
tions of the master equation are also important steps to¬ 
wards the understanding of functional forms of the PDFs. 
Our analysis can be easily extended to three dimensions 
and the extension to other cases is also straightforward, 
e.g. the solutions under shear can be obtained if we apply 
our results for compression and decompression to princi¬ 
pal compressive and tensile directions, respectively jl24| . 
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Appendix A: Discretizations of the PDFs, CPDs, 
and the Chapman-Kolmogorov equation 

In this Appendix, we explain how to calculate the 
PDFs and CPDs from the numerical data. 

We first make a histogram of scaled overlaps before 
(de)compression as hist[£]. Here, we count the number 
of particle pairs, i and j, of which scaled overlaps, 
are in between 


£ - A £/2 < fa < £ + A£/2 , (Al) 

where A£ is a small increment or bin-size for £. Then, 
the PDF of scaled overlaps is defined as 


m) 


hist[£] 

A P A£ ’ 


(A2) 


where N p is the number of Delaunay edges (the total 
number of contacts and virtual contacts) which is given 

by 


Np = hist £] • ( A3 ) 

—oo<^<°° 

From Eqs. (IA2I) and (IA3I) . the normalization condition 
for the PDF, Eq. (jSj). is automatically satisfied as 


E 

— oo<£<°° 


hist[ 5 ] 

N~ P A£ 


A£ = 1 • 


(A4) 
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After (de)compression, we make another histogram of 
new scaled overlaps as hist[£'|£]. Here, we count the num¬ 
ber of particle pairs, i and j, of which scaled overlaps 
after (de)compression, £C, are in between 

- A£'/2 < 4- < £' + A'£/2 , (A5) 

where A£' = A£ is a bin-size for £'. Note that Zij 
were Zij which satisfied Eq. ED before (de) compression. 
The new histogram, hist[£'|£], is connected with hist[£] 
through the following relations: 


hist[£] = 

E 

hist [£'|£] , 

(A6) 


— oo<£'<< 

30 


hist[£'] = 

E~ 

hist [Z'\Z] ■ 

(A7) 


— oo<^<°° 


Then, the joint probability distribution is introduced as 


P(Z';Z) 


N p ■ 


(A8) 


The CPD is defined as t he ra tio of the joint probability 
distribution to the PDF |l!4| . Therefore, the CPD is 
given by using Eqs. (IA2I) and (IA8I) as 


At first, we divide the Chapman-Kolmogorov equation 
© for positive scaled overlaps as 

pOO 

= / wcdg 

J 0 

+ [ w vc (amt;)dz , (bi) 

J —00 

where the first and second terms on the right-hand-side 
represent gains of contacts from contacts (CC) and vir¬ 
tual contacts (VC), respectively. Then, the difference be¬ 
tween the PDFs after and before non-affine deformation 
is given by 

P*+st(0 - P<t>{£,') = 

poo pO 

/ wcc(rmm+ / wwrmm 

J 0 J —00 

-p^a{f o °° Wcc m'n+j ww(£im} 

poo 

= / {Wcc(£'m(o - Wccmwo} dt 

J 0 

+ [ {Wycormto - wcvimwo) (B2) 

J —00 


mm 


p(Z';Z) 

hist[£]A£' 


(A9) 


It is readily found that Eq. (IA9I) satisfies the normaliza¬ 
tion condition for the CPD, Eq. (USD, i-e. 


y- hist[£'|g] , _ 
hist [£] AC ? 


(A10) 


where we used the normalization condition of the CPDs, 
Eq. m- Dividing Eq. (IB2I) by a small increment of 
area fraction, Sip, we find the master equation for positive 
scaled overlaps (£' > 0) as 

Emo = l°° [Tcc(e\z)pm - T CC (z irmo] 

+ f [T V c(z'\z)pm-T CV (z\z')pm)\dz, 

J —00 

(B3) 


where we used Eq. (IA6I) . 

Prom Eqs. (IA2I) and (IA9I) . the Chapman-Kolmogorov 
equation is derived as 


/ OO 

pmwiz'mz 

-00 


y- hist[C] hist^l^] 

oofeoc hiSt ^ A ^' 


1 

N P AZ' 


E hist im 

— °o<^<°° 


hist 

N P AZ' 


(All) 


which yields the discretized PDF after (de) compression, 

p*±h(Z')- 


Appendix B: A derivation of the master equation 

In this Appendix, we derive complete forms of the mas¬ 
ter equation, i.e. Eqs. (l33l) and (13H) . from the Chapman- 
Kolmogorov equation (0 • 


where the transition rates in (CC), (VC), and (CV) are 
defined as 


Tcctm =mm ’ 

84>—t 0 0(p 

(B4) 


(B5) 


(B6) 


respectively. 

Similarly, we divide the Chapman-Kolmogorov equa¬ 
tion for negative scaled overlaps as 

P++6*(t') = f W vv (Z'\Z)P<l>(Z)dZ 

J —00 
p 00 

+ / Wcv(Z'\Z)PmdZ , (B7) 

Jo 

where the first and second terms on the right-hand-side 
represent gains of virtual contacts from virtual contacts 
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(VV) and contacts (CV), respectively. The difference be¬ 
tween the PDFs after and before non-affine deformation 
is given by 

p 4 >+s<i i (£0 - P<t>(€') = 

0 poo 

W vv (e\OP^)dC+ / Wcv(?\OP+{S)dt 

-oo J 0 

-m') { f w vv m')d(i + j™ wvc(*im} 

= [° {Wvvtfizmo-wvvitMPtia}# 

J — oo 
poo 

+ / {Wcv(?\OPt(® - WvcimWO} <%, (B8) 
Jo 

where we used the normalization condition of the CPDs, 


Eq. (l26l) . Dividing Eq. (IB8I) by 5(f>, we find the master 
equation for negative scaled overlaps (£' < 0) as 

h p ^' ] = I- [Tvv ^'^ p ^ ) ~ Tvv ^t'md')} da 

poo 

+ / PcWf m(0 - T VC (dlewd')} dd , 

Jo 

(B9) 

where the transition rate in (VV) is defined as 

Tw(d'\d) = lim Wvv }f 10 ■ (BIO) 

5tj>—>0 0(p 
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